An automated calculation pipeline for differential pair interaction energies with molecular force fields using the Tinker Molecular Modeling Package

An automated pipeline for comprehensive calculation of intermolecular interaction energies based on molecular force-fields using the Tinker molecular modelling package is presented. Starting with non-optimized chemically intuitive monomer structures, the pipeline allows the approximation of global minimum energy monomers and dimers, configuration sampling for various monomer–monomer distances, estimation of coordination numbers by molecular dynamics simulations, and the evaluation of differential pair interaction energies. The latter are used to derive Flory–Huggins parameters and isotropic particle–particle repulsions for Dissipative Particle Dynamics (DPD). The computational results for force fields MM3, MMFF94, OPLS-AA and AMOEBA09 are analyzed with Density Functional Theory (DFT) calculations and DPD simulations for a mixture of the non-ionic polyoxyethylene alkyl ether surfactant C10E4 with water to demonstrate the usefulness of the approach. Scientific Contribution To our knowledge, there is currently no open computational pipeline for differential pair interaction energies at all. This work aims to contribute an (at least academically available, open) approach based on molecular force fields that provides a robust and efficient computational scheme for their automated calculation for small to medium-sized (organic) molecular dimers. The usefulness of the proposed new calculation scheme is demonstrated for the generation of mesoscopic particles with their mutual repulsive interactions.


Introduction
The quantitative description of non-bonding interactions between molecules is fundamental to understanding and designing chemical processes in materials and life sciences [1,2].In contrast to covalent bonding within molecules, non-bonding intermolecular interactions comprise dispersed variations of electromagnetic interactions like dipole/dipole, dipole/induced dipole, induced dipole/ induced dipole (van der Waals) interactions, hydrogen bonding, (partial) charge interactions, π-π, cation/ anion-π or polar π-effects.Each different spatial configuration of two molecules can be assigned a corresponding nonbonding intermolecular interaction energy.Its determination is challenging because the intermolecular interactions are generally small compared to covalent bonding and especially to the total energy of a molecular configuration.
For an accurate quantitative description, a quantum chemical treatment with a suitable model chemistry is commonly advised, e.g., application of Density Functional Theory (DFT) with an appropriate combination of functional and basis set: the complexation energy (i.e., the energy difference between a specific dimer configuration and the two monomer molecules that form it) can then quantify the intermolecular interaction.In particular, Symmetry-Adapted Perturbation Theory (SAPT) allows direct computation of non-bonding intermolecular interactions (i.e., without the need to calculate the total energies of monomers and dimer) and provides a physically meaningful decomposition of its contributing (electrostatics, induction, dispersion, short-range repulsion) terms [2].Recent DFT-SAPT approaches have demonstrated a comparatively fast calculation in combination with remarkable accuracy for small-to medium-sized dimers up to the adenine-thymine base pair [3,4].
When multiple spatial dimer configurations need to be sampled, quantum chemical approaches become increasingly expensive due to their considerable computational complexity, as a single calculation can easily take minutes or even longer.As a purely classical alternative, molecular force fields may be employed instead: they allow intermolecular energy calculations within fractions of a second for a specific spatial dimer configuration, i.e., an acceleration by several orders of magnitude.This comes at the expense of accuracy, as a given force field can easily lead to deceptively erroneous results.
This work aims at providing a robust automated molecular force-field based calculation pipeline for comprehensive estimation of mutual intermolecular energies of a set of small to medium-sized monomer molecules.For a chosen force field, the monomer molecules must be provided with an initial, chemically intuitive spatial geometry with associated atom types (where for small molecules a planar 2D geometry provided by any 2D structure editor seems to be sufficient).Then, each monomer start geometry is globally optimized to its minimum energy conformer.Several monomer-monomer distances with multiple configurations at each distance are sampled to obtain a near-global minimum energy dimer.This dimer is then locally optimized towards its global minimum.The sampled configurations can be averaged by Boltzmann weights to get mean non-bonded dimer interactions at different distances.The calculations can be extended to obtain differential molecule pair interaction energies, which describe the excess intermolecular interaction of two different molecules in comparison to the average interaction of the molecules themselves (see Eqs. 1 and 2 in "Methods" section below), where molecular dynamics (MD) simulations are used to estimate the coordination numbers of neighboring molecules.
The resulting interaction energies may be useful for different purposes.Based on the comprehensive sampling, suitable spatial start configurations can be obtained for more elaborate (e.g., quantum chemical) refinement procedures.Pairwise non-bonded intermolecular interactions can be considered as molecule-pair descriptors for Cheminformatics tasks like molecular similarity estimation.Differential molecule pair interaction energies play an important role in statistical thermodynamics, e.g., for the quantitative estimation of excess quantities that determine the properties of mixtures [5].In particular, they can be related to Flory-Huggins interaction parameters for polymer models (Eq.9), which in turn can be used to describe isotropic particle-particle repulsions for mesoscopic simulation approaches (Eq.10) [6].The latter application is particularly studied in this work.
The backbone of the calculation pipeline is implemented using the Java programming language.All forcefield-based calculations are performed with the Tinker Molecular Modelling package [7].The Tinker package is modularized into stand-alone, task-based executables (marked in italics in "Methods" section below), which fit well into the Java backbone-driven parallelized computational scheme that fully exploits the computational capabilities of multicore workstations.All force fields provided by Tinker can be used for the calculation pipeline.

Methods
The force-field-based intermolecular energy E C ij (r) between two molecules i and j is determined by differ- ent non-covalent-bonding contributions (van der Waals and partial charge interactions, hydrogen bonding etc.) and depends on the intermolecular distance r between the centers of the molecules as well as their relative spatial configurations C .For each specific distance r fix there is a minimum energy configuration C min with E C min ij r fix ≤ E C ij r fix .The global minimum energy dimer is characterized by a distinct distance r min so that r min is the distance between the centers of two molecules i and j when the dimer geometry corresponds to the global energy minimum.If different spatial configurations with a fixed intermolecular distance r fix are averaged, an averaged intermolecular energy for this distance E ij r fix is obtained: The corresponding minimum energy distance r E ,min with �E ij � min = �E ij � r �E�,min ≤ �E ij �(r) does not necessarily coincide with minimum distance r min of the global minimum energy dimer.Among the concrete averaged configurations with a fixed intermolecular distance r fix the configuration C * with the minimal intermolecular energy is denoted A differential pair interaction energy describes the excess intermolecular interaction of molecules i and j in comparison to the average interaction of the molecules themselves.This may be defined with respect to the global minimum energy dimers or corresponding averages (see Eqs. 5 and 6 below for calculation details) A positive (negative) differential pair interaction energy indicates a less (more) favorable intermolecular interaction compared to the individual ones.When differential pair interaction energies are applied to lattice models, all vertices of the lattice have a fixed number of surrounding neighbours (e.g.Z = 4 in two dimensions or Z = 6 in three dimensions).In contrast, continuum models can (and usually do) have different coordination numbers Z ij , where the number Z ij of molecules j surrounding mole- cule i can (and usually does) differ from the number Z ji of molecules i that surround molecule j .Equations 1 and 2 can be extended accordingly to and (1) (3) Thus, the estimation of differential pair interaction energies requires two steps: (I) A calculation scheme to obtain the different (averaged) molecule pair interaction energies E min ij ( E ij min ), and (II) a procedure to estimate the different coordination numbers Z ij .In the follow- ing, the concrete implementation of a corresponding automated calculation pipeline for a selected molecular force field using the Tinker Molecular Modelling package is described.All dimer-related calculations are started with the global minimum energy conformers of the two constituent monomer molecules (see following "Global minimum energy monomers" section), where multiple dimer configurations are analyzed with these global minimum energy monomers (see following "Global minimum energy dimers" section) to obtain an approximate configuration for the global minimum energy dimer.The latter is achieved by optimizing this final approximate dimer configuration without constraints using all atomic degrees of freedom, so that the monomers are no longer confined to their individual global minimum energy conformers.

Global minimum energy monomers
The global minimum energy conformers are derived in advance from conformer search procedures: STARTING from an initially defined chemically intuitive geometry of a monomer molecule, a first geometry improvement is achieved with Tinker optimize using the Optimally Conditioned Variable Metric (OCVM) optimization technique [by default with a root mean square (RMS) gradient of 0.01 kcal/mole/Å] [8].The resulting locally optimized geometry is then transferred to a low-mode conformational search (LMOD) with Tinker scan to find the minimum energy conformer (where by default, the RMS gradient for the LMOD procedure is set to 0.0001 kcal/mole/Å, the energy threshold for local minima is set to 100 kcal/mole, torsion angles are automatically selected, and five eigenvectors are used for the local search) [9].For small molecules with O(10) number of atoms, the detected minimum energy conformers usually coincide with the global minimum energy conformers.

Global minimum energy dimers
To approximate the global minimum energy dimer of a pair of molecules i and j , the centers of the mol- ecules are positioned at different distances ranging (by default) from 3 to 16 Å in steps of 0.5 Å where the initial relative configuration of the two molecules is arbitrary.For each distance, a configuration sampling procedure is performed, which is sketched in Fig. 1.A (unit) sphere around each center is constructed with a number of N sphere evenly spaced points being generated on each sphere using a Fibonacci lattice as an adequate approximation (in comparison to a latitude-longitude lattice the surface points generated by a Fibonacci lattice are more evenly spaced with a smaller axial anisotropy) [10].By rotating the molecules around their centers so that two adjacent spherical surface points and the centers of both molecules lie on a straight axis, N sphere × N sphere configurations are generated for which the corresponding interaction energies are determined by Tinker analyze (with settings to only compute the non-bonding interactions for a significantly accelerated performance).In addition, for each single configuration the second molecule is rotated for a number N rot of angles around the axis with cor- responding interaction energy calculations, so that in total N sphere × N sphere × N rot spatial configurations are sampled for a single fixed distance between the monomer molecules (with each monomer being constrained to its individual global minimum energy conformer).The distance resolution is then refined around the distance with the detected minimum interaction energy dimer by reducing the step size from 0.5 to 0.1 Å, and then again from 0.1 to 0.01 Å (compare Fig. 2), so that a final near minimum energy dimer configuration is evaluated for the latter resolution.This resulting configuration C * is then optimized with Tinker optimize without constraints using all atomic degrees of freedom (i.e. the monomers are no longer confined to their individual Fig. 1 Configuration sampling for the acetic acid dimer: a An initial dimer of two acetic acid molecules is constructed with a distinct distance between the centers of both molecules (with each monomer being constrained to its individual global minimum energy conformer).b Spheres around the centers of the molecules are populated with a number N sphere of evenly distributed (equidistant) points on their surfaces.c The interaction energy is determined for every configuration where two adjacent spherical surface points and the centers of both molecules lie on a straight axis (which is achieved by corresponding rotations of the molecules around their centers).d For each single configuration in c the second molecule is also rotated for a number of angles N rot around the axis with a corresponding interaction energy calculation.e By varying and refining the distance between the molecule centers, the minimum energy dimer is determined and finally optimized to approximate the global minimum energy dimer without constraints using all atomic degrees of freedom.For the acetic acid dimer two perspectives of the finally optimized global minimum energy dimer with E min ij are shown, using the MMFF94 force field [11] global minimum energy conformers) to arrive at the force-field dependent approximation for E min ij of the two molecules i and j , the global minimum energy dimer using the specified force-field.

Configuration sampling
If configuration sampling is desired for a specific distance r fix of the dimer molecules, the (already) obtained inter- action energies of the N sphere × N sphere × N rot configura- tions for this distance may be averaged with Boltzmann weights w C r fix for a defined temperature (with k B , the Boltzmann constant, and T , thermodynamic tempera- ture) using E min ij : Then E ij can be evaluated from the different E ij (r) .Figure 2 shows the quantitative results for the acetic acid dimer using the Merck molecular force field (MMFF94) [11]. (5)

Coordination numbers
The coordination numbers Z ij are estimated by MD simulations.The simulation box construction is based on a pure water box at 298 K to get consistent results.
A water molecule has a van der Waals volume of V H 2 O vdW = 17.35Å 3 , in a pure water box it occupies a vol- ume of V H 2 O box = 30.00Å 3 at 298 K due to its density [12] and molar mass [13].This relation is mapped to other molecules X with so that the edge length a of a cubic simulation box of N molecules X is given by The van der Waals volumes are approximated with the VABCVolume [14] descriptor of the Chemistry Development Kit (CDK) [15,16].A simulation box with a defined number N of molecules j (default is 400) and defined edge length a is created using Tinker xyzedit.Then a single molecule i is added to the box, (7) where Tinker xyzedit automatically removes molecules j to keep the defined density of the simulation box.The (possibly unfavorable) start configuration is optimized with Tinker minimize to avoid atomic contacts that could lead to instabilities.The following MD simulation uses Tinker dynamics with a step size of one femtosecond and an Andersen thermostat [17] for temperature equilibrium (default is 298 K). 10,000 (default) initial steps are used for box equilibration, followed by several hundred thousand steps with data recording (default is 400,000).For each recorded simulation step the number of molecules j surrounding the single molecule i is analyzed.This is done either by a "brute-force" counting approach or, alternatively, by the cell-index method with periodic boundary conditions [18].The "brute-force" approach calculates the distances between each atom of the single molecule i and each atom of all molecules j .Alternatively, the cell-index method only considers the (drastically reduced) distances between each atom of single molecule i and the atoms of solvent molecules j in neigh- bouring cells.The criterion for including a molecule j as a relevant neighbor for the coordination number Z ij is based on the distance between its atoms and those of molecule i .If the distance of the respective atoms is less than or equal to the sum of their van der Waals radii plus an arbitrary "catch radius" (default is 1 Å), the molecules are considered neighbors.In Fig. 3, a snapshot of a simulation step is illustrated.For each simulation step, the number of neighboring molecules j is determined.The average over all recorded steps is used to estimate the coordination number Z ij , see Fig. 3.

Flory-Huggins and mesoscopic repulsion parameters
Differential pair interaction energies may be directly utilized to estimate Flory-Huggins interaction parameters χ ij by to describe polymer solutions [5], with E ij Z being defined in Eq. 4.
For "bridging the gap between atomistic and mesoscopic simulation" (Groot and Warren [6]), the interacting molecules can be identified with the particles of "bottom-up" mesoscopic Dissipative Particle Dynamics (DPD), where (9)  the microscopic Flory-Huggins interaction parameters χ ij can be directly related to mesoscopic isotropic particleparticle repulsions a ij (T ) (expressed in units of k B T , with ρ DPD being the dimensionless DPD density, refer to [6] for details) which determine the conservative DPD forces , soft repulsive DPD force and harmonic bond force on particle i exerted by particle j; a ij , maximum isotropic repulsion between particles i and j ; r ij = r i − r j = r ij r 0 ij ; r 0 ij , unit vector.The numerical factor (3.4965) in Eq. 10 is derived from Eq. 24 in reference [6] where the inverse value (0.286) is given.
In interplay with the dissipative (frictional) forces F D ij and random forces F R ij the conservative forces F C ij guide the trajectories r i (t) of the DPD particles according to Newton's equation of motion: with m i , r i , mass and position vector of particle i ; t , time; F i , total force on particle i exerted by other particles j ; N , total number of particles in simulation; F C ij , F D ij , F R ij , conservative, dissipative, and random force on particle i exerted by particle j.
Thus, the described calculation pipeline may be applied to construct a force-field-based particle set for DPD simulations.

Pipeline code availability
The pipeline code is written in Java and openly available at [19].A dedicated installer executable for the Java pipeline code, which comprises a full Java runtime environment, is available at [20] for the Windows operating system.For Linux operating systems a zip file is available at [20].According to its licensing terms the Tinker executables for optimize, scan, analyse etc. have to be downloaded from its website [21] into a specified pipeline directory: a detailed instruction how to perform this is provided at [22].A set of stand-alone Mathematica notebooks [23] for result visualizations is provided at [24].A test pipeline is available at [19] to ensure proper installation.(10)

Pipeline calculation performance
Calculation of a full single differential pair interaction energy for the force fields MM3, MMFF94 and OPLS-AA with default settings ( N sphere × N sphere × N rot = 144× 144 × 16 = 331,776 dimer configurations for each fixed distance to approximate the intermolecular interaction energies, 10,000 equilibration steps and 400,000 simulation steps for the MD simulations to estimate the coordination number with 400 molecules in the box) takes several hours, where the AMOEBA09 force field requires a multiple.Since the pipeline supports comprehensive calculation parallelization for a set of monomer molecules, a single differential pair interaction energy can be obtained on average in less than an hour on a 64-core AMD Ryzen ™ Threadripper ™ PRO 5995 CPU workstation [25].

DFT calculations for result evaluation
DFT calculations are performed with Gaussian 16 [26] and analyzed with GaussView 6 [27].All molecular geometries are optimized using the dispersion-corrected wB97XD functional [28] with the 6-311++G(d,p) basis set where counterpoise calculations are used to obtain basis set superposition error (BSSE) corrected interaction energies.All Gaussian jobs files used are openly available at [29].

DPD simulations for result evaluation
All DPD simulations of this study are performed with the MFsim simulation system [30] using the Jdpd simulation kernel [31].All constructed particle sets and MFsim simulation jobs are openly available at [32].

Results and discussion
To demonstrate the applicability of the different steps of the calculation pipeline, several small molecules are selected: Water (abbreviated H2O), methane (Me), ethane (Et), methanol (MeOH), dimethyl ether (Me2O) and the acetic acid dimer (HAc).The calculation results for this molecule set are evaluated and compared with alternative approaches and experimental results.All averaged energies and MD simulations refer to a temperature of T = 298K .All intermolecular energy calculations were performed with N sphere × N sphere × N rot = 144 × 144 × 16 = 331,776 dimer configurations for each fixed distance of the molecule centers.

Acetic acid dimer
The acetic acid dimer is stabilized by two hydrogen bonds and has a planar geometry.Calculation results with the MMFF94 force field are shown in Figs. 2 and  4. Figure 2

Mutual dimer calculations
Table 1 presents the mutual dimer calculations: it comprises E ij and E C * ij for the detected minima (compare Fig. 2) as well as the global force-field-based energy minimum E min ij for force fields MMFF94, AMOEBA09, MM3 [34] and OPLS-AA [35,36] with the water models TIP3P [37] and TIP5P [38].As expected, the nonpolar pure alkyl (methane and ethane) dimers exhibit only small interaction energies, the acetic acid dimer with two hydrogen bonds shows the largest interaction, and the polar dimers are in between.There are clear differences between the force fields, with the OPLS-AA (TIP5P) interactions for the alcohol-water dimers being the strongest.On average, MM3 differs most significantly from the other force fields.
For the HAc-HAc and the Me-Me dimer the results for force fields OPLS-AA, AMOEBA09 and MMFF94 were compared with corresponding DFT single point calculations (denoted DFT sp).In addition, the spatial configurations with E min ij were used as start geometries for DFT geometry optimizations (denoted DFT opt), see Table 2.The DFT calculations indicate that the automated pipeline leads to acceptable near minimum energies and corresponding spatial configurations-with individual exceptions: in contrast to the MMFF94 and OPLS-AA force fields, AMOEBA09 results in an eclipsed minimum energy configuration for the Me-Me dimer instead of a staggered one (see Fig. 6, this eclipsed configuration is maintained by the DFT geometry optimization), but this finding has no significant influence on the subsequent investigations due to its low energetic effect (see Table 2).The difference in interaction energies between the DFT opt and the OPLS-AA force field result for the HAc-HAc dimer is most significant.

Coordination numbers Z ij
Table 3 contains the mutual coordination numbers Z ij (at T = 298K ), where a single molecule i is surrounded by 400 molecules j, using the MMFF94, OPLS-AA (with TIP3P and TIP5P water models), and AMOEBA09 force fields.For comparison, static packing [40] results are included which were taken from previous research [41] and obtained with the commercial Blends software Table 1 Force-field based intermolecular interaction energies in kcal/mole for the different dimers (averages are obtained at T = 298K) Dimer MMFF94 MM3 AMOEBA09 OPLS-AA (TIP3P) OPLS-AA (TIP5P)     ) are taken from [29] for acetic acid and [33,39] for methane.For the Me-Me dimer the AMOEBA09 force field yields an eclipsed minimum energy configuration instead of a staggered one (labelled c)  of Materials Studio [42] using the Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies (COMPASS) force field [43].
Figures 7 and 8 show the similarity of trends in coordination number assignment.While the MD-derived coordination numbers are similar for the different force fields used, the results for the COMPASS static packing approach are significantly reduced by about 50%.The spread of the values derived from MD (Fig. 7) is significantly higher than that of the static packing approach.Interestingly, linear mapping of the coordination numbers to the interval [0,1] yields an approximate overlap of the results (Fig. 8).

Repulsion parameters a ij
For the different force fields, particle sets for mesoscopic DPD simulations with the isotropic mutual repulsions a ij for methane (Me), ethane (Et), methanol (MeOH), ethanol (EtOH), dimethyl ether (Me2O), and water (H2O) were constructed.The off-diagonal a ij values of a particle set were linearly scaled with MFsim [30] so that the maximum absolute deviation between the smallest a ij value and the diagonal values a ii = 24.83(for a thermodynamic temperature of 298 K) is 20.For the OPLS-AA force field the water models TIP3P and TIP5P were considered, denoted OPLS-AA (TIP3P) and OPLS-AA (TIP5P).An additional OPLS-AA (TIP5P) particle set, denoted OPLS-AA (TIP5P Z ij = 1 ), was created which is solely based on the minimal averaged intermolecular energy E ij with a fixed coordination number Z ij = 1 for all dimers.With force field MM3, interaction energies can be calculated only, so a combined particle set was created using MM3 for interaction energy calculation and the MMFF94 force field for MD-based coordination number estimation, denoted MM3/MMFF94.The particle set from [41,44], based on the COMPASS force field, is used for comparison.
Figure 9a-f display the different repulsion parameters a ij .The red dashed line indicates the diagonal value of 24.83.A crucial difference between the different particle sets is the water-methanol (H2O-MeOH) repulsion: for the COMPASS, OPLS-AA (TIP5P) and OPLS-AA (TIP5P, Z ij = 1 ) force fields, this repulsion is the smallest one (and thus the base value for off-diagonal repulsion parameter scaling), whereas for force fields OPLS-AA (TIP3P), AMOEBA09, MMFF94, and MM3/MMFF94   force field exhibits an extraordinary difference for the water-methane (H2O-Me) and water-ethane (H2O-Et) repulsions, which is otherwise not visible.A significant difference between the OPLS-AA (TIP5P) and OPLS-AA (TIP5P, Z ij = 1 ) force field is the water-dimethyl ether (H2O-Me2O) repulsion.These obvious differences between the force fields lead to different levels of usefulness for DPD simulation approaches, where even a single difference in dimer interaction can become a decisive factor.

DPD simulations
The created particle sets were used for DPD simulations of mixtures of water with the non-ionic polyoxyethylene alkyl ether surfactant C 10 E 4 , where "C 10 " indicates the number of 10 carbon atoms in the alkyl chain of the lyophobic part, and "E 4 " represents the number of 4 lyophilic ethylene oxide units [44].A stable lamellar L α phase is formed by a C 10 E 4 /water mixture around 298 K with a C 10 E 4 mass fraction of 0.75.The performance of the particle sets for the different force fields may be evaluated by monitoring the emergent formation of C 10 E 4 bilayers from initial random mixing in the simulation box [44].
The DPD simulations are carried out with the settings outlined in [44], using the SPICES 9Me-4Me2O-MeOH fragmentation for C 10 E 4 [45,46], an integration step size of 0.04, and a deactivated periodic boundary in z-direction (to induce bilayer formation in the xy-plane).
For particle set OPLS-AA (TIP5P Z ij = 1 ), a stacked bilayer superstructure emerges at simulation step 62,000, see Fig. 10 and Table 4, which was even below the COMPASS particle set from [44] with 116,000 steps, where the emerged bilayer structure corresponds well to the one reported in [44], see Fig. 11.The OPLS-AA (TIP5P) particle set required more than twice as many simulation steps (132,000) and the OPLS-AA (TIP3P) particle set more than tenfold as many (846,000 steps).The AMOEBA09 and MMFF94 particle sets show a bilayer superstructure formation, but the bilayers do not align parallel to the xy-plane (as induced by the periodic boundary conditions) but parallel to the yz-plane.Using the hybrid MM3/MMFF94 particle set, no bilayer was formed within 1,000,000 simulation steps.For the specified DPD simulation task, the OPLS-AA (TIP5P Z ij = 1 ) particle set can be regarded as the most suitable choice, which is in good agreement with experimental findings.

Conclusions
The outlined automated comprehensive calculation of intermolecular interaction energies based on molecular force-fields shows satisfactory results for small molecule interactions and can even be successfully used to estimate mesoscopic simulation parameters.Special care should always be taken, as individual force fields can lead to erroneous results.Therefore, despite all automation, manual checking of the results is still essential.
The new calculation pipeline can be easily extended to additional force fields (which may require their conversion into the Tinker format).Therefore, calculating differential pair interaction energies will advance with the progress in improving the underlying molecular force fields.Due to the modularized pipeline approach, alternative modeling packages can also be used for certain computational tasks if they can provide the required specific functions.
The robustness of the outlined computational pipeline can be seen as a crucial advance, as the estimation of differential pair interaction energies along alternative computational paths often led to ambiguous results, which in particular prevented the construction of consistent particle-particle repulsions for larger mesoscopic particle sets [47].The construction of mesoscopic sets with dozens of particles (including hundreds or thousands of mutual particle repulsions) is now within practical reach.

Fig. 2
Fig. 2 Acetic acid dimer.Red: E C * ij (r) for the minimum sampled energy configuration, Blue: Averaged interaction energy E ij (r) for temperature T = 298K

Fig. 3
Fig. 3 Averaged coordination number Z ij of a single acetic acid molecule in 400 water molecules starting after 10,000 initial steps (for box equilibration) with a "catch radius" of 1 Å using the MMFF94 force field.The averaged coordination number converges to a value of 17.7.Box graphics: Left: Snapshot of a simulation step of a single acetic acid molecule in 400 water molecules.Right: Magnification of the neighboring water molecules around the single acetic acid molecule.Yellow spheres: atoms of the single acetic acid molecule with their van der Waals radii.Grey spheres: Neighboring water molecules considered for the coordination number determination depicts the minimal energy configuration C * energies E C * ij r fix for each specific distance r fix in red, exhibiting a single minimum at r fix = 4.91 Å with E C * ij (4.91 Å) = − 15.9 kcal/mole.The corresponding averaged intermolecular energies E ij r fix are shown in blue, where the single minimum distance coincides at r fix = 4.91 Å with E ij (4.91 Å) = − 15.4 kcal/mole in this case.The two hydrogen bonds of the minimal energy configuration C * have a length of 1.69 Å and a distance of 2.63 Å between the corresponding donor and acceptor oxygen atoms (see Fig. 4), which is close to the experimental values of 2.68 Å [33].With the final optimization of the near minimal energy configuration C * the global MMFF94 force- field-based minimum energy configuration C min with E min ij = E C min ij (r min ) = −17.6 kcal/mole is obtained, which is 1.7 kcal/mole below the sampled minimal energy C * configuration: the finally optimized dimer keeps the distance of r min = 4.91 Å but shows a planar geometry with a slightly reduced hydrogen bond length

a:
Hydrogen bond length (distance oxygen to oxygen), b: Distance between carbon atoms, c: Eclipsed minimum energy configuration The Me-Me dimer distance denotes the distance between the carbon atoms (b), while the HAc-HAc dimer distance refers to the specific hydrogen bond length (labelled a).Experimental values (denoted Exp.

Fig. 6
Fig. 6 Minimum energy E C * ij configuration for the Me-Me dimer: Left: OPLS-AA force field with staggered configuration.Right: AMOEBA09 force field with eclipsed configuration

Fig. 7 Fig. 8 Fig. 9
Fig. 7 Coordination numbers Z ij , scaled to 1 for each force field.The black line indicates the COMPASS static packing coordination numbers for the homo dimers

Fig. 11
Fig. 11 Me2O and MeOH particle distribution snapshots along the z-axis perpendicular to a single emerged C 10 E 4 bilayer for the OPLS-AA (TIP5P Z ij = 1 ) (solid lines) and COMPASS (dashed lines) particle sets.The highlighted area corresponds to the bilayer width, indicated by the double arrow �E ij � min = �E ji � min but Z ij may be differ- ent from Z ji , i.e.Z ij = Z ji .

Table 2
Force-field based interaction energies E C * ij and E min ij with corresponding DFT interaction energies and configuration measures

Table 3
Coordination numbers Z ij (at T = 298K)